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Abstract: We present the full results of the Wuppertal-Budapest lattice QCD collabora- 
tion on flavor diagonal and non-diagonal quark number susceptibilities with 2-1-1 staggered 
quark flavors, in a temperature range between 125 and 400 MeV. The light and strange 
quark masses are set to their physical values. Lattices with Nt = 6, 8, 10, 12, 16 are 
used. We perform a continuum extrapolation of all observables under study. A Symanzik 
improved gauge and a stout-link improved staggered fermion action is utilized. All results 
are compared to the Hadron Resonance Gas model predictions: good agreement is found 
in the temperature region below the transition. 
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1 Introduction 

The QCD transition, once occured in the Early Universe, is being routinely reproduced 
in the laboratory, in the ultrarelativistic heavy ion collision experiments at CERN SPS, 
RHIC at Brookhaven National Laboratory, ALICE at the LHC and the future FAIR at 
the GSI. The most important known qualitative feature of this transition is its cross-over 
nature at vanishing baryo-chemical potential [1]. A lot of effort has been invested, both 
theoretically and experimentally, in order to find observables which can unambiguously 
signal the transition. As expected in a cross-over, observables follow a smooth behaviour 
over the transition. The characteristic temperature of the transition depends on how one 
defines it. For the renormalized chiral condensate the Wuppertal-Budapest collaboration 
predicted a value around 150 MeV, which was recently confirmed by hotQCD (see for their 
journal publication [2]). 

Correlations and fluctuations of conserved charges have been proposed long ago to 
signal the transition [3, 4]. The idea is that these quantum numbers have a very different 
value in a confined and deconfined system, and measuring them in the laboratory would 
allow to distinguish between the two phases. 

Fluctuations of conserved charges can be obtained as linear combinations of diagonal 
and non-diagonal quark number susceptibilities, which can be calculated on the lattice at 
zero chemical potential [5, 6]. These observables can give us an insight on the nature of 
the matter under study [5, 7]. Diagonal susceptibilities measure the response of the quark 
number density to changes in the chemical potential, and show a rapid rise in the crossover 
region. At high temperatures they are expected to be large, if the quark mass is small 
in comparison to the temperature. At very high temperatures diagonal susceptibilities 
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are expected to approach the ideal gas hmit. On the other hand, in the low-temperature 
phase they are expected to be small since quarks are confined and the only states with 
nonzero quark number have large masses. Agreement with the Hadron Resonance Gas 
(HRG) model predictions is expected in this phase [8]. Non-diagonal susceptibilities give 
us information about the correlation between different flavors. They are supposed to van- 
ish in a non-interacting quark-gluon plasma (QGP). In the approximately self-consistent 
resummation scheme of hard thermal and dense loops Ref. [9] shows nonzero correlations 
between different flavors at large temperatures due to the presence of flavor-mixing dia- 
grams. A quantitative analysis of this observable allows one to draw conclusions about 
the presence of bound states in the QGP [10] . Another observable which was proposed to 
this purpose, and which can be obtained from a combination of diagonal and non-diagonal 
quark number susceptibilities, is the baryon-strangeness correlator [11]. 

Several results exist in the literature about the study of quark number susceptibilities 
on the lattice both for 2 [12] and 2-1-1 [13] quark flavors. However, for the first time 
in this paper the susceptibilities are calculated for physical values of the quark masses 
and a continuum extrapolation is performed not only for strange quark susceptibilities 
[14] but also for the light quark and the non-diagonal ones. We present full results of 
our collaboration for several of these observables, with 2-1-1 staggered quark flavors, in a 
temperature range between 125 and 400 MeV. The light and strange quark masses are 
set to their physical values. Lattices with Nt = 6, 8, 10, 12, 16 are used. Continuum 
extrapolations are performed for all observables under study. We compare our results to the 
predictions of the HRG model with resonances up to 2.5 GeV mass at small temperatures, 
and of the Hard Thermal Loop (HTL) resummation scheme at large temperatures, when 
available. 



2 Observables under study 

The baryon number B, strangeness S and electric charge Q fluctuations can be obtained, at 
vanishing chemical potentials, from the QCD partition function. The relationships between 
the quark chemical potentials and those of the conserved charges are as follows: 

1 2 

I 1 

II , ^ 

= ^fj-B - - I-IS- (2.1) 

Here the small indices u, d and s refer to up, down and strange quark numbers, which, 
too, are conserved charges in QCD. The negativ sign between and /j-s reflects the — 1 
strangeness quantum number of the strange quark. 
Starting from the QCD pressure, 

^ = :^^lnZ(y,r,/iij,/i5,/iQ) (2.2) 
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we can define the moments of charge fluctuations as follows: 
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(2.3) 



In the present paper we will concentrate on the quadratic fluctuations, thus I + m + n = 2. 
In terms of quark numbers {Nx) our observables read^: 
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and on the correlators among different charges or quark flavors: 
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where X and Y are one of u, d and s. Given the relationships between chemical potentials 
(2.1) the diagonal susceptibilities of the conserved charges can be obtained from quark 
number susceptibilities in the following way: 
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If we do not wish to take further derivatives, we can take all three chemical potentials 
(n, d, s) to zero. In this case, for our 2+1 flavor framework nothing distinguishes between 
the u and d derivative: this gives slightly simplified formulae: 
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The baryon-strangeness correlator, which was proposed in Ref. [11] as a diagnostic to 
understand the nature of the degrees of freedom in the QGP, has the following expression 
in terms of quark number susceptibilities: 



Cbs = -3^^ = 1 + (2.8) 



^ For simplicity we inculde the normalization in the definition of X2 a-nd Xii^- In Refs. [14-16] 

we used the notation x$ /T'^ for the same observable. 
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3 Details of the lattice simulations 



3.1 The lattice action 

The lattice action is the same as we used in [15, 16], namely a tree-level Symanzik improved 
gauge, and a stout-improved staggered fermionic action (see Ref. [17] for details). The 
stout-smearing [18] yields an improved discretization of the fermion-gauge vertex and re- 
duces a staggered artefact, the so-called taste violation (analogously to ours, an alternative 
link-smearing scheme, the HISQ action [19] suppresses the taste breaking in a similar way. 
The latter is used by the hotQCD collaboration in their latest studies [2]). Taste symmetry 
breaking is a discretization error which is important mainly in the low temperature phase. 
In the continuum limit the physical spectrum is fully restored. 

In analogy with what we did in Refs. [15, 16], we set the scale at the physical point by 
simulating at T = with physical quark masses [16] and reproducing the kaon and pion 
masses and the kaon decay constant. This gives an uncertainty of about 2% in the scale 
setting. 

For details about the simulation algorithm, renormalization and a discussion on the 
cut-off effects we refer the reader to [16, 20]. 

3.2 Finite temperature ensembles 




Figure 1. The statistics used in this study. The number of trajectories exceeds 10^ for several 
temperatures. Each bar refers to the respective color-coded lattice resolution in a 10 MeV wide 
temperature bin. We analyzed the gauge configurations after every tenth trajectroy with 128 pairs 
of random sources (256 at Nt = 16) with the same physical quark masses that we had in the 
simulation. 

The compact Euclidean spacetime of temperature T and three-volume V is discretized 
on a hypercubic lattice with A'^ and Ng points in the temporal and spatial directions, 
respectively: 

T=^, V = {N,af, (3.1) 
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where a is the lattice spacing. At fixed Nt, the temperature can be set by varying the 
lattice spacing. This implies varying the bare parameters of the lattice action accordingly, 
keeping the pion and kaon (Goldstone) masses at their physical values. In other words, all 
our simulation points lie on the line of constant physics, determined at zero temperature 
in our earlier works [15, 16]. For every given Nt we keep the geometry fixed, such that the 
aspect ratio is ~ 3. 

For the present analysis we use five lattice spacings for each temperature in the transi- 
tion region, corresponding to the temporal resolutions Nt = 6, 8, 10, 12 and 16. We enriched 
our existing set of temperatures since Refs. [14, 21] and four-folded the statistics at Nt = 16. 
We save a configuration every tenth trajectory in a hybrid monte carlo stream with unit 
trajectory length. The statistics used in this analysis is given in Fig. 1, summed in bins of 
10 MeV. 

3.3 Fluctuations from the lattice 

The fluctuations of interest are derivatives of the free energy with respect to the chemical 
potential of a conserved charge. This guarantees the finiteness of the lattice observables, 
thus no renormalization is necessary. 

A derivative of the partition function can be written in terms of S^s, the action with 
all fermionic degrees of freedom already integrated out, as follows: 

dilogZ=^ J VUd^e-^^^ = {Ai) . (3.2) 

When we take further derivatives, the following chain rule applies: 

didj log Z = {AiAj) - {A,) {Aj) + 6ij {B,) . (3.3) 

Here i indicates the variable of the derivative, the chemical potential jii in this case, 
with i = u,d, or s. Ai and Bi are the first and second derivatives of 5eff without the factor 
g-Scff^ Their ensemble averages are calculated with the same weight used for generating 
the configurations. 

In our case Ai and Bi are 

Ai = -^ixMr^M[, (3.4) 
B, = Jtr (MfMri - MiM-'MlMr\) . (3.5) 

with Mi = nrii+Ip is the fermion operator with the bare mass mj, that we also used for 
generating these configurations. M[ and M'-' stand for its first and second derivatives with 
respect to /ij, respectively. The pre-factor ^ is required by the staggered formulation of 
the single flavor trace. These derivatives are mass independent. In the lattice simulation 
as well as in the subsequent analysis, the bare mass is the only parameter that identifies a 
particular type of quark. The B, Q and 5 quantum numbers are provided by Eq. (2.1). 

Ref. [5] describes a stochastic technique for calculating the traces in Eqs. (3.4) and 
(3.5). The traces are rewritten in terms of inner products of random sources. The most 
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expensive part of the present analysis is the calculation of the trace in Eq. (3.4), which 
contains disconnected contributions and appears in almost all susceptibilities as Xud- It 
required 128 pairs of random sources per configuration (256 for Nt = 16). For each pair of 
sources one needs two inversions of the fermion matrix with the light quark mass. 

3.4 Continuum extrapolation 

With five lattice spacings per temperature we are in the position to go beyond the simplest 
form of continuum extrapolation and fit a second order polynomial in a^. Especially at low 
temperatures, such fit is indeed necessary, as the coarser lattices have corrections beyond 
the term. In general, a continuum fit benefits from higher order terms, but this also 
introduces ambiguities, such as whether it is appropriate to keep the coarsest point in the 
continuum extrapolation. The answer to this question is obtained by performing all possible 
extrapolations, and weighting them by the resulting goodness of the fit. Accordingly, we 
varied the number of included lattice spacings and made a linear and a quadratic fit in 
a^. We double the number of such choices by considering extrapolations of the inverse 
fluctuations (1/x) too, and then taking the inverse of the corresponding continuum result. 

There is another source of systematic error: the interpolation ambiguity. The ensem- 
bles were not taken exactly at the same temperatures for different Nt values, and the spline 
fit on the data for a given Nt depends on the node points. We take two choices of the node 
points into the analyses (selecting the original temperature values with either even or odd 
indices). In most cases the two interpolations agree within statistical errors. We incorpo- 
rated the systematic error from this source into the statistical error of the interpolating 
data points prior to the continuum fit. We selected the temperature range for each data 
set such that a consistent interpolation is possible. 

This procedure has (with few exceptions) preferred the full quartic fit over four or five 
points in the transition region, and a suppressed quartic term fit for T > 200 MeV. In most 
cases the reciprocal fit was preferred over the original variant. 

We systematic errors are defined through the central 68.2% of the weighted distribution 
of all analyses, following our collaboration's standard technique c.f. [22, 23]. For simplicity, 
we give our results with the sum of the statistical and the symmetrized systematic errors. 
The continuum bands in the results section correspond to this combined error around a 
central value. In most cases the systematic error dominates. For Xus and Cbs the two 
types of errors are of equal magnitude. The actual smallness of Xud makes the relative 
combined error grow beyond 50%, thus we dropped this observable from our result list. 

We give one example of our interpolation strategy in Fig. 2. We plot the measured X2 
values with statistical errors for two different temperatures. The data points seem to lie 
on a parabola, when plotted as a function of ~ 1 /Nt . We give three fits for each of the 
selected temperatures to indicate the spread of the possible continuum results. 

As discussed in Ref. [21] one can use the tree level improvement program for observ- 
ables. Independently whether one used it or not the results are the same (c.f. Figure 8 of 
Ref. [21]). For simplicity we use in the present paper the direct method and do not ap- 
ply the tree level improvement for our observables when we extrapolate to the continuum 
limit. The improvement factors for the various Nt discretizations (cg = 1.517, cg = 1.283, 
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Figure 2. Examples of our continuum extrapolations. Here we show xu, an observable severely hit 
both by taste breaking and by the cut-off effects in the one-link staggered dispersion relation. The 
data points suggest a quadratic fit in Here we give three possible fits both below and above 

the transition temperature. The solid, dashed and dash-dotted curves represent the fits on the finest 
five, four and three lattices, respectively. The uncentainty related to this choice is incorporated in 
the systematic error (see main text). The statistical errors are much smaller than the size of the 
points, nevertheless the fits provide reasonable values. 

cio = 1.159, ci2 = 1.099, ci6 = 1.054) are merely used here for plotting the raw lattice 
data. 

4 Results 

The first observables we discuss are the diagonal light and strange quark number suscep- 
tibilities: their behavior as functions of the temperature is shown in the two panels of 
Fig. 3. The different symbols correspond to different values of Nt, from 8 to 16. The 
red band is the continuum extrapolation, obtained from the unimproved data, not from 
the improved ones. The continuum extrapolation is performed through a parabolic fit in 
the variable (l/Nt)'^, over five Nt values from 6 to 16. The band shows the spread of the 
results of other possible fits, as discussed in Section 3.4. The comparison between the im- 
proved data and the continuum bands in the figure shows the success of the improvement 
program throughout the entire temperature range. But even the unimproved data could 
be easily fitted for a continuum limit, the combined errors are below 3% in the deconfined 
phase. Both observables show a rapid rise in a certain temperature range, and reach ap- 
proximately 90% of the ideal gas value at large temperatures. However, the temperature 
around which the susceptibilities rise is approximately 15-20 MeV larger for strange quarks 
than for light quarks. In addition, the light quark susceptibility shows a steeper rise with 
temperature, compared to the strange quark one. As expected, they approach each other 
at high temperatures. The effect more evident in Figure 4: in the left panel we show the 
continuum extrapolation of both susceptibilities on the same plot. In the right panel we 
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Figure 3. Left panel: diagonal light quark susceptibility as a function of the temperature. Right 
panel: diagonal strange quark susceptibility as a function of the temperature. In both panels, the 
different symbols correspond to different Nt values. The red band is the continuum extrapolation. 
The black curve is the HRG model prediction for these observables. The dashed line shows the 
ideal gas limit. The light blue band in the left panel is the HTL prediction taken from Ref. [9]. 
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Figure 4. Left panel: direct comparison between the continuum limit of light and strange quark 
susceptibilities. Right panel: ratio X2/X2 ^ function of the temperature. The red band is the 
lattice continuum result. The black, solid curve is the HRG model prediction. The dashed line 
indicates the ideal gas limit. 



show the ratio Xs/Xu- it reaches 1 only around 300 MeV, while for smaller temperatures 
it is < 1. It is worth noticing that all these observables agree with the corresponding HRG 
model predictions for temperatures below the transition. 

The pattern of temperature dependence is strongly related the actual quark mass. The 
difference between the light and strange susceptibilities here with physical masses is more 
pronounced than in earlier works with not so light pions (see E.g. Ref. [13]). 

The non-diagonal us susceptibility measures the degree of correlation between different 
flavors. This observable vanishes in the limit of an ideal, non-interacting QGP. However, 
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Figure 5. Non-diagonal u-s correlator as a function of the temperature. The different symbols 
correspond to different Nt values. The red curve is the continuum extrapolated result. The black 
curve is the HRG model prediction. The dashed line indicates the ideal gas limit for this observable. 

Hard Thermal and Dense Loop framework provides a non-vanishing value for this corre- 
lation also at large temperatures [9]. We show our result in Fig. 5. Xii is non-zero in 
the entire temperature range under study. It has a dip in the crossover region, where the 
correlation between u and s quarks turns out to be maximal. It agrees with the HRG 
model prediction in the hadronic phase. This correlation stays finite and large for a cer- 
tain temperature range above Tc. A quantitative comparison between lattice results and 
predictions for a purely partonic QGP state can give us information about bound states 
survival above Tc [10]. 

Quadratic baryon number, electric charge and isospin fluctuations can be obtained 
from the above partonic susceptibilities through Eqs. (2.7). We show our results for these 
observables in Fig. 6 and in the left panel of Fig. 7. In the low-temperature, hadronic 
phase we have a very good agreement with the HRG model predictions. In the crossover 
region these quantities all show a rapid rise with temperature, in analogy with what already 
observed for the light and strange quark susceptibilities. At large temperature they reach 
approximately 90% of their respective ideal gas values. A comparison between all diagonal 
susceptibilities, rescaled by their corresponding Stefan-Boltzmann limits, is shown in the 
right panel of Fig. 7, from which it is evident that they all show similar features in 
their temperature dependence, even if the temperature at which they rise is larger for the 
strangeness and baryon number susceptibilities. 

The baryon-strangeness correlator Cbs defined in Eq. (2.8) was proposed long ago 
[11] as a diagnostic for strongly interacting matter. It is supposed to be equal to one for 
a non-interacting QGP, while it is temperature-dependent and generally smaller than one 
in a hadronic system. We show our result for this observable in Fig. 8. At the smallest 
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Figure 6. Left panel: isospin susceptibility as a function of the temperature. Right panel: electric 
charge susceptibility as a function of the temperature. In both panels, the different dots correspond 
to different Nt values. The red band is the continuum extrapolation. The black curve is the HRG 
model prediction for these observables. The dashed line shows the ideal gas limit. 
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Figure 7. Left: quadratic fluctuation of baryon number as a function of the temperature. The 
different symbols correspond to different Nt values, the red band is the continuum extrapolation 
and the black, solid curve is the HRG model result. The ideal gas limit is shown by the black, 
dashed line. Right: comparison between all diagonal susceptibilities, rescalcd by the corresponding 
ideal gas limit, as functions of the temperature. 



temperatures it agrees with the HRG model result, and it shows a rapid rise across the 
transition. It reaches the ideal gas limit much faster than the other observables under 
study, yet there is a window of about 100 MeV above Tc, where its value is still smaller 
than one. In analogy with Xiij this observable also gives us information on bound state 
survival above Tc. 

For convenience we tabulate our continuum results in Table 1. 
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Figure 8. Baryon-strangeness correlator as a function of tlie temperature. Tiie different symbols 
correspond to different Nt values, the red band is the continuum extrapolation and the black, solid 
curve is the HRG model result. The ideal gas limit is shown by the black, dashed line. 

5 Conclusions 

In this paper we have presented the continuum results of our collaboration on diagonal and 
non-diagonal quark number susceptibilities, in a system with 2+1 staggered dynamical 
quark flavors with physical masses, in a temperature range between 125 and 400 MeV. The 
continuum extrapolations were based on Nt = 6, 8, 10, 12 and 16 lattices. We calculated 
the systematic errors by varying over the ambiguities of the possible extrapolations. 

All observables consistently show a very good agreement with the HRG model predic- 
tions for temperatures below the transition. 

The diagonal fluctuations have some common features: they all show a rapid rise in 
the crossover region, and reach approximately 90% of the corresponding ideal gas value 
at large temperatures. The rise of both strange quark and baryon number susceptibilities 
is shifted to temperatures about 20 MeV higher than those for light quark, charge and 
isospin susceptibilities. Non-diagonal flavor and charge correlators remain different from 
their ideal gas values for a certain window of temperatures above the transition. This 
pattern encourages further studies to explore the possibility of bound state survival above 
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0.904(15) 


0.900(16) 


0.600(10) 


0.873(15) 


0.294(5) 


-0.0033(8) 


0.9905(22) 


340 


0.908(14) 


0.905(15) 


0.603(9) 


0.882(14) 


0.297(5) 


-0.0028(11) 


0.9920(20) 


360 


0.911(14) 


0.908(14) 


0.605(9) 


0.889(14) 


0.299(5) 


-0.0024(9) 


0.9932(34) 


380 


0.913(15) 


0.911(15) 


0.607(10) 


0.894(16) 


0.300(5) 


-0.0018(10) 


0.9943(39) 


400 


0.915(16) 


0.913(17) 


0.608(11) 


0.899(16) 


0.302(5) 


-0.0012(11) 


0.9953(36) 



Table 1. In this tabic we list the results of our continuum extrapolations. We indicated the sum 
of the statistical and symmetrized systematic errors around the central value. 
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